Analysis of GSE137001

Sys.setenv(`_R_S3_METHOD_REGISTRATION_NOTE_OVERWRITES_` = "false")

suppressPackageStartupMessages({
  library(here)
  library(tidyverse)
})

# Load relevant data
source(here("3_analysis", "data", "get_lfc_dfs.R"))
load(here("3_analysis", "data", "lfc_dfs.RData"))
source(here("3_analysis", "data", "get_enrichr_res.R"))
load(here("3_analysis", "data", "enrichr_res.RData"))

# knit options
options(knitr.duplicate.label = "allow")

Pathway Heatmaps (Enrichr Scores)

# Extract combined scores/p-values from relevant enrichr results
pathways <- c("GO:0009267", "GO:0034198", "GO:0036003", "GO:0036500", "GO:1990928", "GO:0034599", "GO:0034976", "GO:0016236", "GO:0010506")

odds_ratios <- lapply(enrichr_res, function(res) {
  or <- res[["GO_Biological_Process_2021"]] %>%
    filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
    arrange(Term) %>%
    pull(Odds.Ratio)
  
  return(or)
})
odds_ratios <- as.data.frame(odds_ratios)

neglog10_pvals <- lapply(enrichr_res, function(res) {
  pval <- res[["GO_Biological_Process_2021"]] %>%
    filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
    arrange(Term) %>%
    mutate(neglog10 = -log10(Adjusted.P.value)) %>% 
    pull(neglog10)
  
  return(pval)
})
neglog10_pvals <- as.data.frame(neglog10_pvals)

combined_scores <- lapply(enrichr_res, function(res) {
  c_score <- res[["GO_Biological_Process_2021"]] %>%
    filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
    arrange(Term) %>%
    pull(Combined.Score)
  
  return(c_score)
})
combined_scores <- as.data.frame(combined_scores)

row.names(odds_ratios) <- enrichr_res[[1]][["GO_Biological_Process_2021"]] %>%
  filter(str_detect(Term, paste(pathways, collapse="|"))) %>%
  arrange(Term) %>%
  mutate(Term = str_wrap(Term, width = 60)) %>% 
  pull(Term)

row.names(neglog10_pvals) <- row.names(odds_ratios)
row.names(combined_scores) <- row.names(odds_ratios)

Odds Ratio

pheatmap::pheatmap(
  odds_ratios,
  cluster_cols = FALSE,
  scale = "row",
)

-log10(adjusted p-value)

pheatmap::pheatmap(
  neglog10_pvals,
  cluster_cols = FALSE,
  scale = "row",
)

Combined Score

pheatmap::pheatmap(
  combined_scores,
  cluster_cols = FALSE,
  scale = "row",
)

DEGs and Volcano

state_d2_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_d4_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_d6_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_d9_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_d12_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_iPSC_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)

state_ESC_vs_MEF

res <- lfc_dfs[[x]]
padj_cutoff <- 0.01
lfc_cutoff <- 0.5

EnhancedVolcano::EnhancedVolcano(
  res,
  lab = res$SYMBOL,
  x = "log2FoldChange",
  y = "padj",
  pCutoff = padj_cutoff,
  FCcutoff = lfc_cutoff,
  subtitle = sprintf(
    "Adjusted p-value cutoff = %.2f\n|Log2FoldChange| cutoff = %.1f",
    padj_cutoff, lfc_cutoff
  )
)

res %>%
  select(SYMBOL, log2FoldChange, padj) %>%
  filter(padj <= padj_cutoff & abs(log2FoldChange) >= lfc_cutoff) %>%
  arrange(padj) %>% 
  DT::datatable() %>%
  DT::formatSignif(columns = 2:3, digits = 3)